// QUIRKSIM
// Copyright (c) 2009 Gregory Andrew Robinson (garobinson@ucdavis.edu)
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.

#ifndef QUIRKSIM_H_
#define QUIRKSIM_H_

#include <string>
#include <vector>
#include "CLHEP/Vector/ThreeVector.h"

// A macro to disallow the evil copy constructor and operator= functions
// This should be used in the private: declarations for a class
#define DISALLOW_COPY_AND_ASSIGN(TypeName)    \
  TypeName(const TypeName&);                  \
  void operator=(const TypeName&)

// A macro to disallow all the implicit constructors, namely the
// default constructor, copy constructor and operator= functions.
//
// This should be used in the private: declarations for a class
// that wants to prevent anyone from instantiating it. This is
// especially useful for classes containing only static methods.
#define DISALLOW_IMPLICIT_CONSTRUCTORS(TypeName) \
  TypeName();                                    \
  DISALLOW_COPY_AND_ASSIGN(TypeName)

using CLHEP::Hep3Vector;

const double PI = 4.0*atan(1.0);

class ParamPoint;
class QuirkSystem;

typedef ParamPoint Point_T;
typedef std::vector<Point_T*> String_T;

// Print a Hep3Vector's value with 20-digit precision to stdout.
void Print(const Hep3Vector& v);

template <class ForwardIterator>
void STLDeleteContainerPointers(ForwardIterator begin,
                                ForwardIterator end) {
  while (begin != end) {
    ForwardIterator temp = begin;
    ++begin;
    delete *temp;
  }
}

// Represents a point's derivative with respect to time.
struct XPDeriv {
  // Derivative of spatial coordinate with respect to time.
  Hep3Vector dx_dt;

  // Derivative of momentum with respect to time.
  Hep3Vector dp_dt;
};

// A point in space with position and momentum.
class XPPoint {
 public:
  // x: Original position.
  // p: Original momentum.
  XPPoint(const Hep3Vector& x, const Hep3Vector& p);

  virtual ~XPPoint();

  // Increment position and momentum by val*factor.
  void Add(const XPDeriv& val, double factor) {
    x_ += val.dx_dt * factor;
    p_ += val.dp_dt * factor;
  }

  Hep3Vector* mutable_x() { return &x_; }
  Hep3Vector* mutable_p() { return &p_; }
  const Hep3Vector& x() const { return x_; }
  const Hep3Vector& p() const { return p_; }

 private:
  // Position.
  Hep3Vector x_;
  // Momentum.
  Hep3Vector p_;

  DISALLOW_COPY_AND_ASSIGN(XPPoint);
};

// A point along the connecting string.
class ParamPoint : public XPPoint {
 public:
  // Just pass along these to XPPoint(Hep3Vector&, Hep3Vector&);
  ParamPoint(const Hep3Vector& x, const Hep3Vector& p);

  virtual ~ParamPoint() = 0;

  // Return a pointer to a new deep copy.
  virtual ParamPoint* Clone(const QuirkSystem& sys) const = 0;

  // Computes momentum and position derivatives for each point in the system.
  // params:
  //   previous: derivatives from previous intermediate step.
  //   dt: timestep.
  //   out: where to write the result.
  virtual void ComputeDerivs(const QuirkSystem& in_state,
                             String_T::const_iterator it,
                             XPDeriv* out) const = 0;

  virtual Hep3Vector Pi() const = 0;
};

// Base class for vector-valued forces as functions of position and velocity.
class ForceInterface {
 public:
  virtual Hep3Vector Force() const = 0;
  virtual ~ForceInterface() {};

 protected:
  ForceInterface() {}
};

class Quirk : public ParamPoint {
 public:
  Quirk(const Hep3Vector& x, const Hep3Vector& p, double mass,
        const QuirkSystem& sys, const bool lambda);
  virtual ~Quirk();
  virtual ParamPoint* Clone(const QuirkSystem& to_sys) const;
  virtual void ComputeDerivs(const QuirkSystem& in_state,
                             String_T::const_iterator it,
                             XPDeriv* out) const;
  Hep3Vector DxDt() const;
  Hep3Vector OutwardUhat() const;
  Hep3Vector VPar() const;
  Hep3Vector VPerp() const;
  Hep3Vector DpDt() const;
  Hep3Vector Pi() const { return VPerp() * sqrt(1.0 / (1 - VPerp().mag2())); }
  void RegisterForce(ForceInterface* f);
  double mass() const { return mass_; }
  const QuirkSystem& sys() const { return sys_; }
  bool lambda() const { return lambda_; }
 private:
  std::vector<ForceInterface*> forces_;
  double mass_;
  const QuirkSystem& sys_;
  bool lambda_;
  DISALLOW_COPY_AND_ASSIGN(Quirk);
};


class FixedEnd : public ParamPoint {
 public:
  FixedEnd(const Hep3Vector& x);
  virtual ParamPoint* Clone(const QuirkSystem& sys) const;
  virtual void ComputeDerivs(const QuirkSystem& in_state,
                             String_T::const_iterator it,
                             XPDeriv* out) const;
  Hep3Vector Pi() const { return Hep3Vector(0, 0, 0); }
};

class StringPoint : public ParamPoint {
 public:
  StringPoint(const Hep3Vector& x, const Hep3Vector& p, const QuirkSystem& sys);
  virtual ~StringPoint();
  virtual ParamPoint* Clone(const QuirkSystem& to_sys) const;
  virtual void ComputeDerivs(const QuirkSystem& in_state,
                             String_T::const_iterator it,
                             XPDeriv* out) const;
  Hep3Vector DxDt() const;
  Hep3Vector DpDt(String_T::const_iterator it) const;
  Hep3Vector Pi() const { return p(); }
  const QuirkSystem& sys() const { return sys_; }
 private:
  const QuirkSystem& sys_;
  DISALLOW_COPY_AND_ASSIGN(StringPoint);
};

class QuirkSystem {
 public:
  QuirkSystem(double sigma, double point_separation);
  ~QuirkSystem();

  // Return pointer to a deep copy.
  QuirkSystem* Clone() const;

  // One fourth-order Runge Kutta step.
  // params:
  //   dt: timestep.
  void RK4Step(const double dt);

  int ConstrainP() {
    std::vector<ParamPoint*>::iterator it;
    for (it = points_->begin();
         it < points_->end();
         it++) {
      *(*it)->mutable_p() = (*it)->p().perpPart(Uhat(it));
    }
    return 0;
  }

  // Returns length of string and recomputes if necessary.
  double Length() const;

  // Returns the length along string between the min-th and max-th points.
  double DeltaL(String_T::const_iterator min,
                String_T::const_iterator max) const;
  // Computes derivatives for each point in system.
  // params:
  //   previous: derivatives from previous intermediate step.
  //   dt: timestep.
  //   out: where to write the result.
  static void ComputeDerivs(const QuirkSystem& in_state,
                            std::vector<XPDeriv>* out);


  // Calculates change in position with respect to length of string.
  Hep3Vector Uhat(String_T::const_iterator it) const;

  const std::vector<ParamPoint*>& points() const { return *points_; }
  std::vector<ParamPoint*>* mutable_points() { return points_; }
  unsigned int NumPoints() { return points_->size(); }

  double sigma() const { return sigma_; }

  friend class Quirk;

 private:
  // Adds to position and momentum, each multiplied by a scaling factor.
  // params:
  //   add: values to be added
  //   factor: multiply each x, p by this before adding, such as a timestep.
  void Add(const std::vector<XPDeriv>& add, double factor);

  // Reparametrizes string such that no span is longer than 2 * separation and
  // no string is shorter than separation / 2.
  void Reparam(double separation);

  // Parameterization of points along string.
  // First and last objects must be quirks.
  String_T* points_;

  double length_;
  double separation_;
  const double sigma_;

  DISALLOW_COPY_AND_ASSIGN(QuirkSystem);
};

// Represents the gauge string interaction for a single quirk in a pair.
class QuirkForce : public ForceInterface {
 public:
  QuirkForce(const Quirk& p);
  Hep3Vector Force() const;
 private:
  const Quirk& quirk_;
  DISALLOW_COPY_AND_ASSIGN(QuirkForce);
};

#endif  // QUIRKSIM_H_
